First, use this code to download the data from the USGS site.

library(dataRetrieval)
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.4     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(lubridate)
library(tidymodels)

── Attaching packages ────────────────────────────────────── tidymodels 1.3.0 ──
✔ broom        1.0.7     ✔ rsample      1.2.1
✔ dials        1.4.0     ✔ tune         1.3.0
✔ infer        1.0.7     ✔ workflows    1.2.0
✔ modeldata    1.4.0     ✔ workflowsets 1.1.0
✔ parsnip      1.3.1     ✔ yardstick    1.3.2
✔ recipes      1.2.0     
── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
✖ scales::discard() masks purrr::discard()
✖ dplyr::filter()   masks stats::filter()
✖ recipes::fixed()  masks stringr::fixed()
✖ dplyr::lag()      masks stats::lag()
✖ yardstick::spec() masks readr::spec()
✖ recipes::step()   masks stats::step()

library(tsibble)

Registered S3 method overwritten by 'tsibble':
  method               from 
  as_tibble.grouped_df dplyr

Attaching package: 'tsibble'

The following object is masked from 'package:lubridate':

    interval

The following objects are masked from 'package:base':

    intersect, setdiff, union

library(plotly)

Attaching package: 'plotly'

The following object is masked from 'package:ggplot2':

    last_plot

The following object is masked from 'package:stats':

    filter

The following object is masked from 'package:graphics':

    layout

library(feasts) #autoplot visualizations

Loading required package: fabletools

Attaching package: 'fabletools'

The following object is masked from 'package:yardstick':

    accuracy

The following object is masked from 'package:parsnip':

    null_model

The following objects are masked from 'package:infer':

    generate, hypothesize

library(fable)

#If you prefer to keep the date as a date object, you can use lubridate::year() and lubridate::month() to extract the year and month, and then combine them:



poudre_flow <- readNWISdv(siteNumber = "06752260",    # Download data from USGS for site 06752260
                          parameterCd = "00060",      # Parameter code 00060 = discharge in cfs
                          startDate = "2013-01-01",   # Set the start date
                          endDate = "2023-12-31") |>  # Set the end date
  renameNWISColumns() |>                              # Rename columns to standard names (e.g., "Flow", "Date")
  mutate(Date = as.Date(Date, format = "%Y-%m")) |>  # Create a Year-Month column as "YYYY-MM"
  group_by(Date) |>                         # Group by the new Year-Month column
  summarise(Flow = mean(Flow, na.rm = TRUE))           # Calculate the average daily flow for each month

GET:https://waterservices.usgs.gov/nwis/dv/?site=06752260&format=waterml%2C1.1&ParameterCd=00060&StatCd=00003&startDT=2013-01-01&endDT=2023-12-31

Convert to tsibble: Use as_tsibble() to convert the data.frame into a tsibble object. This will allow you to use the feast functions for time series analysis.

# Assuming you already have 'poudre_flow' with the "Date_YearMonth" column as shown earlier
poudre_flow_tsibble <- poudre_flow %>%
  as_tsibble()  # Convert to tsibble using Date as the time index

Using `Date` as index variable.

Plotting the time series: Use ggplot to plot the time series data. Animate this plot with plotly

p_plot<- poudre_flow_tsibble %>%
  autoplot() +
  geom_line(color= "purple") +
  labs(
    title= "Interactive time series plot of the Poudre Rive flow",
    x= "Date",
    y= "Flow")

Plot variable not specified, automatically selected `.vars = Flow`

ggplotly(p_plot)

3. Subseries: Use gg_subseries to visualize the seasonal patterns in the data. This will help you identify any trends or seasonal cycles in the streamflow data. Describe what you see in the plot. How are “seasons” defined in this plot? What do you think the “subseries” represent?

#make month column to predict seasons + group by them
poudre_flow_tsibble <- poudre_flow_tsibble %>%
  mutate(month = month(Date)) %>%
  group_by(month)

gg_subseries(poudre_flow_tsibble, Flow) +
  labs(title= "Seasonal Patterns in Poudre River Stream Flow",
       y= "Stream Flow",
       x= "month",
       color= month) +
  theme(
    text = element_text(size = 10),             # Reduce font size
    axis.text.x = element_text(angle = 45, hjust = 1),  # Rotate x-axis text if needed
    strip.text = element_text(size = 9)          # Resize month labels (facets)
  ) +
  theme_minimal()

The series shows a default of 12 months in the cycle, however it’s very disorganized. If it was beginning from January (1) than the trend shows that in the winter the flow is lower, then begins to pick up in spring months, before falling in summer and autumn months. Some years have better flow in their months. The subseries represents the months within the data for each year.

Decompose: Use the model(STL(…)) pattern to decompose the time series data into its components: trend, seasonality, and residuals. Chose a window that you feel is most appropriate to this data. Describe what you see in the plot. How do the components change over time? What do you think the trend and seasonal components represent?

# STL decomposition model doesn't assume constant seasonal effect like decompose()

# Convert to a regular time series object for STL (monthly data)
flow_ts <- ts(poudre_flow_tsibble$Flow, start = c(2013, 1), frequency = 12)

# Apply STL decomposition
stl_decomposed <- stl(flow_ts, s.window = "periodic") %>%
  plot()

This plot shows the seasonal trends from 2013 and predicts them up until about 2350. They follow typical oscillations that occur on a yearly basis, except they slowly start having lower maximum flows. The trend could represent the effects of climate change on the river flows over time, and the seasonal component is driven by the typical natural biological processes as well as climate change (some seasons may start to get slightly shorter).